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In self-organized critical (SOC) systems avalanche size distributions follow power-laws. 
Power-laws have also been observed for neural activity, and so it has been proposed 
that SOC underlies brain organization as well. Surprisingly, for spiking activity in vivo, 
evidence for SOC is still lacking. Therefore, we analyzed highly parallel spike recordings 
from awake rats and monkeys, anesthetized cats, and also local field potentials from 
humans. We compared these to spiking activity from two established critical models: the 
Bak-Tang-Wiesenfeld model, and a stochastic branching model. We found fundamental 
differences between the neural and the model activity. These differences could be 
overcome for both models through a combination of three modifications: (1 ) subsampling, 
(2) increasing the input to the model (this way eliminating the separation of time scales, 
which is fundamental to SOC and its avalanche definition), and (3) making the model 
slightly sub-critical. The match between the neural activity and the modified models held 
not only for the classical avalanche size distributions and estimated branching parameters, 
but also for two novel measures (mean avalanche size, and frequency of single spikes), and 
for the dependence of all these measures on the temporal bin size. Our results suggest 
that neural activity in vivo shows a melange of avalanches, and not temporally separated 
ones, and that their global activity propagation can be approximated by the principle that 
one spike on average triggers a little less than one spike in the next step. This implies 
that neural activity does not reflect a SOC state but a slightly sub-critical regime without 
a separation of time scales. Potential advantages of this regime may be faster information 
processing, and a safety margin from super-criticality, which has been linked to epilepsy. 



Keywords: self-organized criticality, human intracranial recordings, spike train analysis, highly parallel recordings, 
spiking neural networks, multiunit activity, cortex, monkeys 



INTRODUCTION 

Avalanches, earthquakes, and forest fires are all cascades of activ- 
ity in otherwise quiescent systems (Gutenberg and Richter, 1944; 
Bak et al, 1987; Drossel and Schwabl, 1992; Frette et al, 1996; 



Measures, variables, and abbreviations: a, connection strength or synaptic 
strength; /3, scaling exponent (DFA); a, branching parameter; <T*, estimated 
branching parameter; t, critical exponent of the avalanche size distribution; bs, 
bin size; DFA, detrended fluctuation analysis; f(s), avalanche size distribution; f(s — 
1 , bs) , frequency of avalanches of size s — 1 and their dependence on the bin size; h, 
rate of input spikes, also called drive (Hz); <s>, mean avalanche size; <IEI>, aver- 
age inter event interval; <IEI> — l/R; JV, number of sampled (model) neurons; 
r, rate per unit (Hz); R, population rate (Hz); STS, separation of time scales. 



Dickman et al., 2000). Most of the time, the size of these cas- 
cades, or avalanches, is small, but sometimes avalanches are large 
enough to span the entire system. The size s of an avalanche is 
the number of units activated during a cascade, and interestingly, 
the distribution f(s) of avalanche sizes in the systems mentioned 
above precisely follows a power law: 

m ~ s- r (i) 

where x is the critical exponent. Critical exponents determine 
the macroscopic behavior of a system, and indicate the system's 
universality class (Wilson, 1975). 
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Power law distributions are characteristic for second-order 
phase transitions, where the system is in a "critical" state. If the 
system evolves to reach a critical state without fine-tuning of 
control parameters, the system is termed self- organized critical 
(SOC) (Bak et al, 1987; Jensen, 1998; Nagler et al, 1999; Beggs 
and Plenz, 2003; Frigg, 2003; Beggs and Timme, 2012; Pruessner, 
2012). 

SOC models show avalanches or cascades of activity across 
their units, which may arise from simple local interactions (Bak 
et al, 1987; Drossel and Schwabl, 1992; Olami et al, 1992). These 
avalanches can include all units in the system. However, most 
avalanches are small or intermediate in size. Note that avalanches 
of size one, i.e., only one unit is active and no further activity 
is triggered, have the highest chance of occurring (see Equation 
1). Overall, avalanches are not characterized by an average size, 
i.e., the size distribution is scale-free, and only the true size of the 
system restricts the avalanche size range. 

In nervous systems, scale-free properties have been observed 
in local field potentials (LFP), electro- and magnetoencephalo- 
graphic (EEC, MEG) activity, and BOLD signals (Linkenkaer- 
Hansen et al, 2001; Beggs and Plenz, 2003; Petermann et al., 
2009; Hahn et al, 2010; Ribeiro et al, 2010; Tetzlaff et al, 2010; 
Friedman et al, 2012; Poil et al, 2012; Tagliazucchi et al., 2012; 
Priesemann et al, 2013; Shriki et al., 2013). They have been found 
in different preparations, ranging from cultures to in vivo prepa- 
rations, and across different species and phyla: leeches, rats, cats, 
monkeys, and humans (Linkenkaer-Hansen et al, 2001; Beggs 
and Plenz, 2003; Mazzoni et al, 2007; Pasquale et al, 2008; 
Petermann et al., 2009; Priesemann et al, 2009, 2013; Hahn et al., 
2010; Ribeiro et al, 2010; Tetzlaff et al, 2010; Friedman et al, 
2012; Poil et al, 2012; Tagliazucchi et al, 2012; Shriki et al, 2013). 
The prevailing hypothesis is that scale-free neural activity arises 
from SOC behavior (Linkenkaer-Hansen et al, 2001; Beggs and 
Plenz, 2003; Mazzoni et al., 2007; Beggs, 2008; Pasquale et al, 
2008; Petermann et al, 2009; Shew et al, 2009; Hahn et al, 2010; 
Ribeiro et al, 2010; Tetzlaff et al, 2010; Friedman et al., 2012; 
Poil et al, 2012; Tagliazucchi et al, 2012; Gal and Marom, 2013; 
Shriki et al, 2013). However, there are also studies that reported 
deviations from scale-free activity: Neural activity was shown to 
exhibit sub-critical and super-critical behavior during develop- 
ment in vitro (Pasquale et al., 2008; Tetzlaff et al., 2010; Friedman 
et al, 2012); and there are also studies in which in vivo neural 
activity appeared as sub-critical (Bedard et al., 2006; Priesemann 
et al., 2013). Thus, healthy brains seem to be capable of organizing 
themselves into a range of states that are not necessarily SOC. 

Nevertheless, because neural activity from coarse scale mea- 
sures (e.g., population spikes, LFP, MEG, BOLD) often do 
show power law scaling, the same was expected for more 
basic constituents of neural activity, namely the spiking activity. 
Surprisingly, however, spike avalanches often deviated from 
power law scaling (Bedard et al, 2006; Pasquale et al., 2008; 
Hahn et al, 2010; Tetzlaff et al, 2010). In fact, to the best of our 
knowledge, there is not a single study that demonstrated power 
laws for spikes in awake animals. The deviations from power law 
scaling in previous studies were attributed either to sub- or super- 
critical states (Pasquale et al, 2008; Tetzlaff et al., 2010), or to 
subsampling effects (Ribeiro et al., 2010). Subsampling refers to 



the technical constraint that only a fraction of all neurons in a 
given area can be measured. Subsampling can impede the obser- 
vation of power law distributions in SOC models (Priesemann 
et al, 2009, 2013; Ribeiro et al, 2010; Girardi-Schappo et al, 
2013) and hence a critical system can be misinterpreted as sub- or 
super-critical (Priesemann et al., 2009). Therefore, subsampling 
effects need to be taken into account when interpreting spike 
avalanches. 

An important property of SOC systems, which is potentially 
absent in neural activity, is the separation of time scales (STS) 
(Bak et al, 1987; Drossel and Schwabl, 1992; Clar et al, 1996; 
Dickman et al., 2000; Pruessner, 2012; Hartley et al, 2013) 
whereby pauses between avalanches last much longer than the 
avalanches proper. For example, forest fires last for a much shorter 
time than it takes to regrow the forest. Similarly, earthquakes 
are much more rapid than the time it takes to build shear stress 
through plate tectonics (Drossel and Schwabl, 1992; Clar et al., 
1996, 1999; Baiesi and Paczuski, 2004). Likewise, in the classical 
sandpile model, scale-free avalanche distributions are observed 
only if the grains are dropped at a low enough rate (Vespignani 
and Zapperi, 1997, 1998). This low rate of external input, called 
drive, is a necessary condition for the long pauses and hence for 
SOC (Bak et al, 1987; Drossel and Schwabl, 1992; Clar et al., 1996; 
Dickman et al, 2000; Pruessner, 2012; Hartley et al, 2013). 

Neither the neural activity we analyzed here, nor that from 
previous studies of neural avalanches showed STS: There were 
no long pauses in the neural activity which could be seen as 
natural separations between avalanches. Without such pauses, 
unambiguous detection of the beginning and the end of an indi- 
vidual avalanche is not possible. Hence, the method of temporal 
binning had been introduced as a workaround (Beggs and Plenz, 
2003) (Figure 1). Here, the choice of the bin size determines what 
is considered to be a pause between avalanches. Consequently, 
avalanche sizes necessarily change with the choice of the bin size 
(see e.g., Beggs and Plenz, 2003; Priesemann et al., 2009, 2013; 
Hahn et al., 2010). This implies that also the avalanche size distri- 
butions and, more importantly, power law exponents change with 
the choice of bin size (Beggs and Plenz, 2003; Priesemann et al., 
2013). This is in marked contrast to fully sampled SOC systems, 
in which the power law exponents do not change under temporal 
binning as a result of STS. These differences have to be considered 
when comparing neural activity in vivo to that of classical SOC 
models. 

As indicated above, in classical SOC systems each avalanche 
is separated from the next one by a long pause. In contrast, in 
driven SOC systems, i.e., SOC systems without STS, avalanches 
can meet, merge, intermingle, and split up: They form a melange. 
As we demonstrate in this paper, neural activity indeed resembles 
such a melange of avalanches instead of well-separated ones. 

To investigate the differences between in vivo and model activ- 
ity, we analyzed spike avalanches recorded in awake rats and mon- 
keys, anesthetized cats, and LFP avalanches recorded in humans, 
and compared these in vivo avalanches to avalanches from an 
established SOC model (Bak-Tang-Wiesenfeld model) (Bak et al., 
1987; Dunkelmann and Radons, 1994; Priesemann et al, 2009, 
2013), and to those from a stochastic branching model (Harris, 
1963; Haldeman and Beggs, 2005). 
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FIGURE 1 | Definition of avalanches sizes, branching parameter a", 
and their change with bin size. (A) To define avalanches, temporal 
binning (boxes) is applied to a sequence of spikes (red dots and 
diamonds). Empty bins are marked in blue. An avalanche is an ensemble 
of spikes in a sequence of non-empty bins. Its size s is the total 
number of spikes, as indicated above the bins. The branching parameter 
or* is the ratio between the number of spikes in one bin, divided by the 
number of spikes in the previous bin, as indicated below the bins. If the 
previous bin was empty, o\ is "not defined" (nd). The estimated 
branching parameter or* for an experiment is the average over all a*, a" 



= <af>. (B) When increasing the bin size, the observed avalanches can 
become larger, since pauses "disappear'.' The branching parameter a* 
also changes with the bin size. (C) Under subsampling, only a fraction of 
the units are recorded (red dots), while others are missed (gray). Thereby 
subsampling can split a single avalanche into several parts. (A-C) In the 
model, spikes are either triggered externally by some drive (red 
diamonds), or they are evoked by presynaptic activity (red dots). If a 
second avalanche is triggered while the first one is still active [last 
avalanche in (A)], then the two avalanches cannot be told apart and are 
evaluated as if they were a single one. 



RESULTS 

As a widely held belief states that mammalian nervous sys- 
tems operate in a SOC state, we first briefly recapitulate 
the theoretically expected avalanche statistics in this state by 
example of a SOC model and a critical stochastic branch- 
ing model. We then show that all of the analyzed neural 
avalanches in vivo showed clear deviations from the expected 
statistics. 

The remainder of the results then demonstrates how two sim- 
ple and neurophysiologically well-motivated conceptual changes 
in the models can serve to align model and in vivo activity with 
respect to a large set of measured quantities. 

DIFFERENCES BETWEEN NEURAL DYNAMICS IN VIVO AND SOC 

The first example model is a simple neural network model, which 
is known to have SOC properties (Bak et al., 1987). Furthermore, 
this SOC model has been shown to match LFP avalanches in mon- 
keys and humans (Priesemann et al., 2009, 2013). In our study, 
the model consisted of 2500 non-leaky integrate-and-fire neurons 
arranged as a 50 by 50 grid with nearest neighbor connections 
of synaptic strength a = 1 (see Methods). In this model, spikes 



are either evoked by activity from presynaptic neurons, or by a 
random external input to a neuron. This input is termed drive 
and has a rate h. For h — > 0 and a = 1, this model obeys local 
energy conservation (Bonachela et al, 2010), and is equivalent 
to the well-known SOC Bak-Tang-Wiesenfeld model (Bak et al., 
1987). h — >■ 0 is necessary for a model to be SOC (Vespignani 
and Zapperi, 1997, 1998; Dickman et al., 2000), because it guar- 
antees the obligatory STS. h — >■ 0 is implemented by applying 
external input only when there is otherwise no activity in the 
model. The input triggers an avalanche, i.e., a cascade of events. 
The size s of an avalanche is defined as the total number of spikes 
evoked by a single input spike. This model is known to show a 
power law for f(s) with slope r ~ 1 (Figure 2A), and a cutoff at 
s « 1000 (BaketaL, 1987). This cutoff reflects the finite size of 
the model (Bak et al, 1988; Kadanoff et al, 1989; Ktitarev et al, 
2000). 

To later demonstrate that our conclusions are not specific to 
the SOC model above, we simulated a second model, namely 
a stochastic branching model (see Methods) (Harris, 1963; 
Haldeman and Beggs, 2005). Like the SOC model, it was sim- 
ulated with 2500 neurons, but in contrast to the SOC model, 
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FIGURE 2 | Avalanche size distributions f(s) changed with the bin size for 
the in vivo spike trains (D-F), and for the subsampled models (B,C,H,I). 

(A) f(s) of the SOC model under full sampling did not depend on the bin size. 

(B) Under subsampling (W = 100 neurons), f(s) of the same SOC model 
changed with small bin sizes only. (C) In the driven model (h > 0) f(s) changed 
for all bin sizes, h was chosen such that the population rate Ft of the 100 



sampled model neurons matched R of the experiments (R 320 Hz). (D-F) 
f(s) recorded in the hippocampus (awake rat), the visual cortex (anesthetized 
cat), and the prefrontal cortex (awake monkey). (G-l) shows the same as 
(A-C), but for a critical branching model instead of the SOC model. Dashed 
lines indicate potential power law slopes to guide the eye. All f(s) are 
logarithmically binned and f(s) is in absolute counts. 



the k = 4 postsynaptic neurons were chosen randomly at each 
step. Activity propagated stochastically, i.e., an active neuron acti- 
vated each of its k postsynaptic neurons with probability p = 
a/k. Like the SOC model, this model is critical for a = 1, and 
sub- (super-) critical for a < 1 (a > 1). The critical stochastic 
branching model with STS also showed a power law distribution 
for f(s), but with a different critical exponent (r = 1.5, 
Figure 2G). 

The results for the stochastic branching model and the SOC 
model were qualitatively the same for all measures used below. 
The similarity also held when the models were modified analo- 
gously. Therefore, in the following, we mainly report results for 
the SOC model. 



Our critical models were contrasted with highly parallel 
recordings from awake rats (hippocampus), awake monkeys (pre- 
frontal cortex), and from an anesthetized cat (visual cortex 
area 18). The avalanche distributions f(s) from these in vivo 
spike recordings were all very similar, but clearly differed from 
those obtained from the fully sampled critical models (compare 
Figures 2D-F with A,G). In particular, the in vivof(s) neither fol- 
lowed a power law, in contrast to what is expected for a SOC 
system, nor an exponential distribution, as would be expected 
for independent Poissonian activity (Figures SI and S2 show the 
in vivo f(s) for each experiment in double-logarithmic and log- 
linear scales, respectively). Quantitatively, the f(s) were best fit in 
16 out of 17 experiments by a lognormal distribution 
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(fa(s)-M) 2 

f(s) ~ e 

with parameters /x = 0.89 ± 0.25 and variance a 2 = 1.2 ± 0.1, 
given a bin size of 1 average inter event interval (<IEI>) (see 
Clauset et al., 2007; Priesemann et al., 2013 for details). Based 
on these parameters the maximum of f(s) was at s = 0.87 ± 
0.38 (mean ± SD), which means that f(s) was monotonically 
decreasing. Two alternative distributions, namely stretched expo- 
nentials and power laws with cutoff, also provided reasonable 
fits, with likelihoods ~1% worse than the one for the lognormal 
distribution. 

Interestingly, all in vivo avalanche distributions were similar 
despite changes in the population rate _R by a factor of 50 (from 
37 to 1560 Hz) across the 17 experiments (Figures SI, S2). 

Note that some of the f(s) of the rat experiments could also be 
approximated by a power law, but at most for one selected bin 
size (Figure S3 A). When slightly changing the bin size, the f(s) 
clearly deviated from power law scaling (Figure S3B). This is in 
stark contrast to the behavior expected for SOC systems. 

A second striking difference between critical models and 
in vivo activity was that the in vivo f(s) changed with the bin size 
across a range from 0.5 to 128 ms. The reason for the bin size 
dependence was that in vivo recordings showed pauses of variable 
length between the spikes, while SOC activity showed only the 
long pauses between avalanches, which are due to STS. In order 
to introduce pauses of variable length into the model avalanches, 
one can apply subsampling and drop STS (see next two sections). 

SUBSAMPLING INTRODUCES PAUSES AT SHORT TIME SCALES 

Subsampling refers to the problem that we are far from being 
able to sample all spikes from all neurons, even for a single brain 
area (Figure 1C). Thus, for a careful comparison between in vivo 
recordings and models, the activity from the models should be 
subsampled in the same manner as in the experiments. Because 
in each experiment around 100 neurons were recorded in paral- 
lel, in the model we constrained the sampling also to N = 100 
randomly chosen neurons out of the 2500. We fixed the subsam- 
pling by the number of neurons, and not the fraction, because 
running these critical models with millions of neurons is beyond 
our computational capacities, and because the qualitative results 
did not change in larger models, i.e., when decreasing the fraction 
(see below). 

When applying subsampling, the model avalanche size distri- 
bution f(s) changed with bin size (Figures 2B,H). A change in 
bin size affected /(sj, because subsampling introduces apparent 
pauses in a single avalanche (Figure 1C). These apparent pauses 
were relatively short compared to the duration of an avalanche, 
and compared to the pauses between avalanches on the full model 
(by definition of STS). Therefore, when subsampling, f(s) changed 
only with small bin sizes but stopped to change its shape with 
larger ones (Figures 2B,H). 

These results also held when using a larger model and sampling 
the same number of neurons, i.e., a smaller fraction of neurons. In 
this case, the distance and hence the traveling time of avalanches 
between sampled neurons became larger and longer, and the inter 
spike intervals became unrealistically long. Nonetheless, at large 
bin size, a similar fraction of small avalanches was observed (due 



to STS). As a consequence, f(s) also stopped changing like in 
smaller models, and never became as flat as the in vivof(s). Hence, 
the behavior of a larger model was the same as that of smaller 
ones, but on a longer time scale. 

Subsampling the SOC model did not only introduce a depen- 
dence oi f(s) on the bin size, it also affected the cutoff of/(sJ. 
Thereby, the absolute value of the cutoff became more similar for 
the model and the in vivof(s) (Figures 2B,H). 

In sum, acknowledging subsampling effects in the model 
allowed for a better match between the model and the in vivo 
activity, but only for small bin sizes up to a few milliseconds. For 
larger bin sizes, the in vivof(s) continued to become flatter, while 
the model f(s) stopped to change their shape. This indicated that a 
modification to the model dynamics itself was necessary to match 
in vivo activity. 

AN INCREASE IN DRIVE RATE h CREATES A MELANGE OF AVALANCHES 

We hypothesized that in vivo and SOC activity differed because 
SOC models have STS (Vespignani and Zapperi, 1997, 1998; 
Dickman et al., 2000), which is necessarily absent in vivo. STS 
can be eliminated from the models by increasing the drive rate 
h. We increased h in such a way that the model population 
rate R matched the in vivo population rate under subsampling 
(h = 0.02 Hz, and R = 320 Hz; single neuron rate r in the model: 
r = R/N = 3.2 Hz). In this driven SOC model, the avalanches 
were no longer separated by long pauses (Figure 3B). Instead, at 
any point in time, avalanches could start, meet, intermingle, split 
into branches, or die out (Figures 1, 3B). In such a melange of 
avalanches, single avalanches can no longer be tracked. 

The melange of avalanches in the driven model hardly showed 
any pauses when all neurons were sampled (Figure 3B). However, 
under subsampling, pauses were more frequent. Thus, sub- 
sampling allowed for an extraction of apparent avalanches by 
applying temporal binning (Figure 1). Note that these appar- 
ent avalanches do not correspond to the avalanches observed 
in classical SOC models in which avalanches are separated by 
long pauses, and are thereby defined unambiguously. However, 
the apparent avalanches from the driven models are conceptu- 
ally the same as those extracted from in vivo recordings because 
avalanches in both cases are extracted with the same method. 

As expected for the driven, subsampled SOC model, f(s) 
changed with all bin sizes (Figures 2C,I), and thereby resem- 
bled the in vivo f(s) much better than the original SOC model 
(Figure 2). 

DRIVEN CRITICAL AND DRIVEN SUB-CRITICAL STATES 

In the following, we address the question whether subsampling 
and the elimination of STS is sufficient to match the model 
activity with the in vivo activity, or whether it is necessary to 
introduce in addition deviations from criticality. 

To tune models away from criticality, we made use of the fact 
that SOC and branching models are only critical in the con- 
servative limit (a = 1) (Harris, 1963; Bonachela and Munoz, 
2009; Bonachela et al., 2010). Hence, by introducing dissipation 
(a < 1) these models can be made sub-critical. In fact, the model 
dynamics showed a smooth transition from the "driven SOC" 
state (a = 1) to pure Poisson activity (a = 0) (Figures 3, 4) with 
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FIGURE 3 | The population spike rate of the (modified) SOC model 
depended on the connection strength a and the rate of input spikes h 
(drive), h and a were balanced such that the rate of each unit was r = 5 Hz, 
except for (A), where a = 1 and h— > 0 (SOC model). In (A), the broken 
axes indicate that the pauses between subsequent avalanches are much 
longer than the avalanche proper (separation of time scales). (B) a = 1 , 
h = 0.02 Hz, r = 5 Hz (driven SOC). (C) a = 0.95, h = 0.5 Hz, r = 5 Hz 
(driven sub-critical). (D) a = 0, h = r = 5 Hz (Poisson). In (A-D), the 
population rate time course is indicated in black; the scale bar indicates the 
firing rate per neuron. Black dots show the spike raster from 100 randomly 
chosen units; the blue background denotes pauses, i.e., none of the 2500 
neurons spiked. Note the absence of pauses in (CD). 



decreasing a. In principle, a decrease in a also decreased the fir- 
ing rate r of each unit. To still maintain a constant firing rate r, 
a concomitant increase in the drive rate h was applied. In this 
way, the model could make the transition from driven SOC to 
Poissonian activity without a change in r (Figure 4, black line). 
Given a fixed r, a decrease in a decreased the variability of the 
models population rate (Figure 3). 

To understand which network dynamics between driven criti- 
cal and Poissonian accounted best for the in vivo spike avalanches, 
we identified those measures in the model which depended most 
sensitively on a under subsampling: a had a prominent effect 
on the avalanche size distribution f(s), in particular how f(s) 
depended on the bin size. We quantified this below using the 
following avalanche measures: the mean avalanche size (<s>), 
the frequency of avalanches of size s = 1 (/(s = 1)), and the esti- 
mated branching parameter a*. The way in which these measures 
changed with the bin size depended sensitively on a. In addition, 
we estimated the scaling exponent fi of the "detrended fluctua- 
tion analysis" (DFA) (Peng et al, 1994, 1995; Kantelhardt et al., 
2002). (Note that the scaling exponent (/j) is often denoted as a 
in the literature). The results of these analyses are presented in 
detail below, and compared one by one to the in vivo results. 




drive h (Hz) 



FIGURE 4 | In the model, the spike rate r of a unit depended on the 
synaptic strength a and the rate of input spikes (/?). With increasing /) or 
a, the rate of each unit increased. The black line indicates the parameter 
combination of a and h, for which r = 5 Hz. 



THE MEAN AVALANCHES SIZE 

The mean avalanche size (<5>) from the subsampled model fol- 
lowed a power law with increasing bin size for a = 1 (driven 
SOC), and followed an exponential for a = 0 (Poissonian activ- 
ity) (Figure 5A). For intermediate values of a, the relation 
changed gradually. 

For the experiments, the observed <s> at a given bin size 
depended strongly on the population spike rate R that varied con- 
siderably between experiments (R ranged from 37 Hz to 1.5 kHz). 
To diminish the impact of R, we used a normalized bin size, i.e., a 
bin size in units of average inter-event-intervals ( 1 <IEI> = l/R). 
Using the normalized bin size, the <s> of all experiments over- 
lapped (Figure 5A, gray lines). However, the <s> did not follow a 
power law with changing bin size in vivo, in contrast to the driven 
critical model. In fact, the in vivo <s> was best matched by the 
<s> of the driven, sub-critical models (a ~ 0.99). Thus, com- 
paring the in vivo and model <s> indicated that spike avalanches 
resembled a driven sub-critical regime more closely than a driven 
SOC state. 

THE FREQUENCY OF AVALANCHES OF SIZE ONE 

The frequency of avalanches of size s = 1, /(s = 1, bs) quanti- 
fies how/(sJ decayed with the bin size (bs) at s = 1, i.e., how the 
intercept oif(s) with the y-axis in Figure 2 changed. f(s) at s = 1 
was equally spaced from bin size 1 to 32 ms for the driven critical 
models under subsampling (Figures 2C,I) which is remarkable as 
it corresponds to a power law behavior of f(s =1, bs) for the 
driven SOC model (black line in Figure 5B; note that the x-axis 
here is in <IEI>, and 1 <IEI> = 2 ms in the model). For the 
sub-critical models (a < l),f(s= 1, bs) decayed more steeply 
than a power law. For the Poissonian case (a = 0), it followed an 
exponential. In this respect, /(s =1, bs) and <5> showed similar 
behaviors with a. 

f(s = 1, bs) is a promising new measure to assess criticality 
under subsampling, because in contrast to many other measures, 
its behavior did not change with the subsampling strategy: For the 
driven SOC model, it showed power law scaling independently 
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FIGURE 5 | Two new avalanche measures. (A) The mean avalanche size and 
(B) the frequency of avalanches with size s = 1 , f(s = 7, fas), changed with the 
bin size (fas) in the model (colored) and in the experiments (gray). The colored 



lines show f(s = 7, fas) for the model with varying synaptic strength a. In the 
model, the drive rate fa was adjusted such that each neuron spiked with r f» 5 Hz. 
In (B), f(s = 1, fas) was normalized such that f(s = 7, fas = 7 </£/>) = 1. 




FIGURE 6 | The frequency of single events f(s= 1, bs). Decreased with 
the bin size (fas) as a power law, independently of the subsampling set 
in the driven SOC model (a — 1 , r — 5 Hz). The subsampling set is indicated 
in the right part of the figure. It was chosen as follows: blue f(s = 1 , fas): 
sampling 64 random units; green f{s= 1, fas): sampling 100 random units 
(both, blue, and green units together); red and turquoise: sampling 8x8 
units arranged in a grid with distance 1, and distance 5, respectively; pink 
and beige: sampling 4x4 units with distance 1, and 5, respectively. 



of the number and spatial arrangement of the sampled units 
(Figure 6). However, the slope of the power law did change 
due to the model's next-neighbor topology: With smaller dis- 
tances between sampled sites, the power laws became flatter 
(red and pink traces in Figure 6). For the stochastic branch- 
ing model, the same results held, but the power law slopes did 
not change under subsampling, owing to the model's random 
topology. 

The in vivof(s = 1 , bs) did not follow a power law (Figure 5B, 
gray lines), and for most cases did not follow an exponential 
dependency either (Figure 5B). The best approximation for the 
in vivo f(s = 1 , bs) was the driven, slightly sub-critical model 
(a ~ 0.99). This is in agreement with the results for <s>. 

The precise value of a necessary to achieve the best match 
between model and experiments potentially depended on a 
number of factors (e.g., finite size effects). However, the main 



result that <s> and/(s =1, bs) observed in vivo followed nei- 
ther a power law nor an exponential distribution excludes both, 
critical and Poissonian states of operation. 

THE BRANCHING PARAMETER a 

A widely used measure to estimate whether the in vivo avalanches 
reflected a driven SOC brain state is the branching parame- 
ter a*, which has been used in many past studies about neural 
avalanches to test whether the brain was SOC (Beggs and Plenz, 
2003; Beggs, 2007; Plenz and Thiagarajan, 2007; Priesemann et al., 
2009, 2013; Shew et al, 2009; Klaus et al, 2011; Shriki et al, 
2013). The analysis of a* was initially inspired by the theory of 
branching processes (Harris, 1963), in which a = 1 guarantees 
that a branching process is critical. Note, however, that estimating 
a* from data may yield misleading results, because a* depends 
on various factors such as the bin size (Beggs and Plenz, 2003; 
Priesemann et al., 2013), the subsampling geometry (Priesemann 
et al, 2009), and STS (i.e., h — > 0 vs. h > 0). We next show how 
a* depended on these factors in our models, and then use these 
results to estimate whether the in vivo avalanches might reflect a 
SOC state. 

For the modified SOC model, we expect that a equals a. For 
the second model we used, i.e. the stochastic branching model, 
we know by definition of the model that a equals a. However, 
when estimating a* in this model by applying temporal binning 
to the model activity, finding the expected a* = a was the excep- 
tion, not the rule (Figure S4; results were very similar to the ones 
for the SOC model in Figure 7). In addition, er* changed with 
the bin size, although the model parameter a proper is obviously 
bin size independent (Figures 7, S4). Although the estimated a* 
failed to approximate the true a ,a* may still be a viable approach 
to compare model and in vivo activity in the following. Since the 
results for both models were basically the same, we again focus on 
the results for the modified SOC model. 

With STS, a* always approached zero for large bin sizes inde- 
pendently of model state and subsampling approach (dashed lines 
in Figures 7A,B, S4). For intermediate bin sizes and under sub- 
sampling, a* varied widely, a* tended to be smaller for smaller 
a, but the absolute value of a* apparently cannot serve as an 



Frontiers in Systems Neuroscience 



www.frontiersin.org 



June 2014 | Volume 8 | Article 108 | 7 



Priesemann et al. 



Spike avalanches in vivo 




bin size (ms) bin size (ms) bin size (<IEI>) 



FIGURE 7 | The estimated branching parameter a* changed with bin 
size. (A,B) In the model, a* depended on the synaptic strength a and the bin 
size. For the driven model, the spike rate was fixed to r = 5 Hz (full lines), 
while for the model with separation of time scales the drive was infinitesimal 
small (h -> 0; dashed lines). For h -> 0 and a = 1 , the model is SOC (black 



dashed lines). (A) Results for the fully sampled model. (B) Results for 
subsampling N = 100 neurons from the model. (C) a* for the spiking activity 
recorded in monkeys, cats, and rats varied with the bin size, but was very 
similar across species and experiments. It was well approximated by the 
driven model with a = 0.98 (green line). 
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FIGURE 8 | The exponent /? of the DFA. Depended on the synaptic 
strength a in the model (diamonds), and was affected by subsampling 
(black: fully sampled model; green: subsampled model). For the 
experiments, fi (gray circles) and the respective mean values (gray bars) 
ranged between 0.55 and 0.9. 



indicator for the state of the system (Figures 7A,B)- Thus, under 
most analysis conditions, the estimated a* did not show the 
intended result (a* = a). Note that in theory, a* should not 
change at all with the bin size. 

Without STS (full lines in Figures 7A,B, S4), a* was <1 for 
small bin sizes, > 1 for intermediate bin sizes, and approximated 
unity for large bin sizes - independently of the state of the model. 
This shows that the widely held assumption that an estimated 
a* > 1 (a* < 1) corresponds to a super-critical (sub-critical) 
state of the system is likely incorrect, especially for the ubiquitous 
scenario of subsampling. 

Although the expected a* = 1 is neither unique to critical 
systems, nor indicative of criticality, a* and its dependence on 
the bin size still reflect the intrinsic dynamics of the system. 
Therefore, comparing a* between in vivo and model activ- 
ity may still help to indicate the state of the system. Note 
that to estimate the in vivo a* we used the normalized bin 
size (in <IEI>) to account for the different population rates 
R in the experiments, er* was very similar across all experi- 
ments (Figure 7C) despite a 50-fold difference in R. This indi- 
cates once again that neural avalanches in vivo hardly dif- 
fer across mammalian species (from rats to monkeys), across 
brain structures (from hippocampus to prefrontal cortex), and 
across cognitive states (from anesthetized to awake behaving 
animals). 

Given the complex dependence of a* on the bin size, how can 
a* be used to estimate the precise state of the neural network? 
First, for all in vivo avalanches, a* approximated unity for large 
bin size (Figure 7C). However, this simply indicates that spiking 
activity in vivo lacks STS. Second, the maximum of a* under sub- 
sampling may be an indicator of the state. The maximum of a* 
increased with increasing a. For a = 1, a* showed a maximum 
of ~3 at bs ~100ms. [The same values held for the stochastic 
branching model (Figure S4)]. For the experiments, the maxi- 
mum value of a* was only around 1.4. Overall, the best match 
for the in vivo a* was achieved by the driven, slightly sub-critical 
models (a ~ 0.98). This result is in line with the previous results 
for/(s = 1, bs) and <s>. 



THE SCALING EXPONENT /? 

In DFA, the scaling exponent ft quantifies the memory decay 
in a time series. /3 = 0.5 indicates that a time series has no 
memory (uncorrelated); /3 ~ 1 indicates 1/f (pink) noise; and 
P ~ 1.5 Brownian noise. We estimated /j for the population rate 
time series of the model (r = 5 Hz), and for each experiment. 
As expected, under full sampling the model with a = 1 showed 
P ~ 1 (Figure 8, black diamonds); with decreasing a, /3 decreased 
as well; and for a = 0 (Poisson), we found fi ~ 0.5. Qualitatively, 
the same results held under subsampling, but fi tended to be 
underestimated (Figure 8, green diamonds). 

The in vivo activity showed neither /3 = 1 nor f) = 0.5, but ft 
ranged between 0.55 and 0.9. These fi values correspond to those 
of the sub-critical, driven model with 0.98 < a < 0.999. 

All the above measures indicated that driven, slightly sub- 
critical models provided the best match to in vivo spike 
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FIGURE 9 | Avalanche size distributions f(s) for in vivo spikes and for 
spikes from the driven, sub-critical models. (A) f(s) of one awake 
monkey. Dots indicate the raw f(s), while lines are the f(s) with 
logarithmic binning. (B) f(s) for the driven, sub-critical models with 
a = 0.99, and r = 5 Hz; model 1 denotes the modified SOC model (full 



lines), and model 2 the stochastic branching model (dashed lines). (C) 
f(s) of all in in vivo spike recordings (rat, cat, monkey), together with the 
f(s) of the driven, subcritical models (same as in B). All bin sizes were in 
average inter event intervals (<IEI>), and f(s) were normalized such that 
f(s= 1, bs= 7):= 1. 



avalanches. Most of these measures were derived from the 
avalanche size distribution, and hence we expect a good match 
between the in vivo f(s), and the f(s) of the driven models with 
a ~ 0.99. Indeed, given a normalized bin size, both sub-critical 
models fitted the in vivof(s) well (Figure 9). The small differences 
for large s (s > 100) may potentially be overcome by applying 
a more realistic drive instead of uncorrelated Poissonian drive, 
for example one that reflects the statistics of neural activity (as 
lined out here), or the statistics of our environment (Field, 1987; 
Van der Schaaf and van Hateren, 1996; Simoncelli and Olshausen, 
2001; Sinz etal, 2009). 

LFP AVALANCHES IN HUMANS 

Approximate power law distributions have been reported for 
coarse measures of neural activity, such as population spikes, 
LFP, EEG, MEG, and BOLD activity (Linkenkaer-Hansen et al., 
2001; Beggs and Plenz, 2003; Petermann et al, 2009; Hahn et al., 
2010; Ribeiro et al, 2010; Tetzlaff et al, 2010; Friedman et al., 
2012; Poil et al, 2012; Tagliazucchi et al., 2012; Priesemann 
et al., 2013; Shriki et al., 2013). In the following, we show 
that also LFP recordings in humans indicate a driven, slightly 
subcritical regime, despite their approximate power law scaling 
of /•'>-,•• 

LFPs were recorded using intracranial depth electrodes from 
five human subjects. Each subject had between 44 and 63 record- 
ing contacts implanted. From these recordings, we extracted 
avalanches of enhanced activity (see Methods and Priesemann 
et al., 2013). The LFP f(s) closely followed a power law 
(Figure 10A), and the slope of the power law decreased with 
increasing bin sizes. This is in contrast to SOC systems in which 
the slope does not change with temporal binning (Figures 2A,G), 
and indicates that LFP avalanches, like the spike avalanches, lack 
clear STS. 

In general, the LFP f(s) showed a better approximation to 
power law scaling than any of the spike avalanche distribu- 
tions (Figures 2, 10). Despite an approximate power law scaling 
for f(s), all the other measures we used here [i.e., <s>, f(s = 



1, bs), a*, and /j] indicated a sub-critical regime: The <s> 
and the f(s= 1, bs) both deviated from power law scaling 
(Figure 10B); the branching parameter did not show a pro- 
nounced peak ( Figure 11); and the scaling exponent /j of the DFA 
was smaller than unity (mean(/3) = 0.6; Figure 7). This is in line 
with our previous study on the same data (Priesemann et al., 
2013), and with our results for spiking activity. In sum, despite 
approximate power-law scaling in/(s), all the other measures indi- 
cated a driven, slightly sub-critical regime on the level of LFP 
activity. 

DISCUSSION 

This study challenges the hypothesis that mammalian brains 
operate in a SOC state, as has been repeatedly suggested 
(Linkenkaer-Hansen et al., 2001; Beggs and Plenz, 2003; 
Haldeman and Beggs, 2005; Levina et al, 2007a; Hsu et al., 
2008; Pasquale et al., 2008; Stewart and Plenz, 2008; Petermann 
et al, 2009; Priesemann et al, 2009; Shew et al., 2009; Hahn 
et al., 2010; Ribeiro et al., 2010; Tetzlaff et al, 2010; Poil et al, 
2012; Tagliazucchi et al, 2012; Shriki et al, 2013). Despite 
these claims, evidence for SOC was found lacking for spik- 
ing data, which are generally considered an important and 
reliable marker of neural activity. To test the SOC hypothe- 
sis, we therefore analyzed in vivo spiking activity from three 
mammalian species and local field potential recordings from the 
human brain using established measures of criticality, and also 
novel ones that are robust to common shortcomings of exper- 
imental data, such as subsampling. We particularly focused on 
systematic changes of these measures with the choice of the 
bin size. 

Spike avalanches from rats, cats, and monkeys, and LFP 
avalanches from humans showed deviations from the behavior 
expected for SOC, thereby contradicting the SOC hypothesis. 
To reproduce the in vivo results and provide potential explana- 
tions for their deviations from SOC, we modified the models 
capable of critical behavior. We found a close match between 
in vivo and model behavior ( 1 ) if those models were subsampled, 
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FIGURE 10 | (A) The size distribution f(s) of LFP avalanches from intracranial 
depth electrodes in humans followed power laws. The slope of the power 
laws changed with the bin size (see legend). The bin size was changed over 
a 1000-fold range, from sampling resolution (400 Hz, i.e., 2.5 ms) to 
"gluing" everything together at bs « 2500 ms. The bin size closest to one 
inter event interval is marked in purple (bs= 80 ms, see Methods). (B) 
Neither the mean avalanche size (<s>), nor the frequency of avalanches of 
size s = 1 , f(s = 1, bs), showed a power law. Each line represents the 
results for one recording session (<s> in black, f{s = 7, bs) in gray). 



and (2) if the STS - a fundamental property of SOC sys- 
tems - was eliminated, and (3) if the models were tuned to a 
sub-critical regime. As these results generalized over two very 
different models, we interpret results from the in vivo record- 
ings here as evidence that mammalian nervous systems operate 
in a driven, sub-critical regime. This regime, albeit not critical, 
was, however, remarkably similar across species and experimental 
conditions. 

UNIVERSAL BEHAVIOR OF SPIKE AVALANCHE DISTRIBUTIONS 
ACROSS RECORDING AREAS. VIGILANCE STATES AND SPECIES 

The observed avalanche size distributions f(s) were similar across 
species and recording areas (hippocampus in rats, visual cortex 
in cats, prefrontal cortex in monkeys). A similar universality of 
f(s) across recording areas has been reported by Ribeiro and col- 
leagues (hippocampus, somatosensory cortex, and visual cortex 
in rats) (Ribeiro et al., 2010). Thus, avalanche activity seems to 
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be independent of the function and the precise anatomy of an 
area. This might either indicate that avalanches are not a sen- 
sitive measure of neural dynamics, or that activity propagation 
must follow principles that are independent of the specific role 
that a brain area plays in information processing. The first argu- 
ment is not likely applicable, since avalanches change under data 
shuffling and they sensitively reflect the correlation structure in 
the data (e.g., Figure 1 in Priesemann et al, 2013). The second 
argument might indeed hold. Hence, the challenge is to identify 
the principle that gives rise to these apparently universal spike 
avalanche distributions. This principle may in fact be very simple. 
As discussed below, our modified SOC model, as well as a simple 
branching model, suggests that on average one spike gives rise to 
a little less than one subsequent spike, and that quiescence in the 
population activity is prevented by "input spikes" which trigger 
avalanches at a low rate. This principle differs from SOC, where 
one spike on average gives rise to exactly one subsequent spike, 
and the rate of input spikes approaches zero (STS). As a conse- 
quence, SOC activity shows only one avalanche at a time, while 
the driven, slightly sub-critical regime shows instead a melange of 
avalanches. 

EMPIRICAL AVALANCHE DISTRIBUTIONS RULE OUT THE CRITICAL 
AND THE P0ISS0N STATES 

Let us first summarize the conclusions that can be drawn from 
the analyses of the in vivo spike avalanches alone, without refer- 
ring to modeling. For f(s), neither was the power law scal- 
ing found, that is characteristic for SOC, nor did the novel 
measures (f(s= 1, bs), <s>) support the hypothesis of criti- 
cal behavior. Thus, the hypothesis that spike avalanches show 
signs of SOC can be ruled out. In addition, we can rule out 
the hypothesis of largely independent Poissonian behavior of 
the spiking units (that is often used in models), because in 
this case the avalanche distributions should have shown expo- 
nential behavior, which was not observed. We therefore con- 
clude that spiking activity is neither (self-organized) critical nor 
Poissonian. 
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LIMITATIONS OF THE MODELS AND MEASURES 

The SOC model used here was admittedly simple - it comprised 
neither inhibitory connections nor leakage in the neurons; synap- 
tic connections had a homogeneous nearest-neighbor topology 
and were all of identical strength a. We chose this model because 
the basic variant (a = 1, h —> 0; i.e., the Bak-Tang-Wiesenfeld 
model; Bak et al., 1987) is extensively studied in the context of 
SOC (De Menech et al, 1998; Jensen, 1998; Vespignani et al., 
1998; Dickman et al, 2000; Dhar, 2006; Pruessner, 2012). The 
second model we used was a stochastic branching model (Harris, 
1963; Haldeman and Beggs, 2005). It was set up to be compa- 
rable to the SOC model, but had a random topology, and the 
activity propagated stochastically with p = a/k. In this model, 
the number of connections k hardly affected the results (see also 
Haldeman and Beggs, 2005). 

For both models, the avalanche dynamics was qualitatively 
similar. Hence, the model results were not specific to the topology 
(local vs. random), the number of connections k, and the pre- 
cise spike propagation mechanisms (deterministic vs. stochastic). 
In contrast, implementing leaky model neurons may hinder SOC 
altogether (Bonachela and Munoz, 2009; Bonachela et al, 2010). 
This in itself is an argument against the hypothesis that neural 
activity is SOC, but it could still be "quasi-critical" (Bonachela 
and Munoz, 2009; Bonachela et al, 2010). However, our results 
indicate sub-criticality. 

We note that the power law scaling observed for the novel mea- 
sures (f(s= 1, bs), <s>) in the critical models has not been 
derived analytically yet. However, in both critical models the 
novel measures showed power law scaling despite the different 
topology and the different spike propagation rules, and hence 
we expect this behavior to be characteristic for critical dynamics. 
Still, for now these measures can only be used as tools to compare 
model and in vivo dynamics, and not for determining scaling laws. 

ON THE PLAUSIBILITY OF EXTERNAL DRIVE 

Spike and LFP avalanches recorded in rats, cats, and primates 
were best matched by a driven sub-critical model. The drive in 
the model consisted of input spikes, i.e., of spikes not caused by 
presynaptic spikes from within the model. Given their impor- 
tance for a successful match between in vivo and model activity, 
we may ask what the in vivo counterpart of the input spikes in 
the models could be. In vivo, such input spikes can be provided 
by at least three sources — by sensory input elicited by stimuli 
in the outside world, from brain structures other than the one 
under consideration, or by internal activation which presumably 
occurs spontaneously. Such spontaneous activity can for exam- 
ple be generated by pacemaker cell activity (Selverston, 2008; 
Longtin, 2013), or vesicle fusion at a presynaptic terminal with- 
out a preceding spike (Fredj and Burrone, 2009). With all these 
known input sources in vivo, it came as no surprise that the model 
required input spikes (i.e., drive) to be able to match in vivo 
activity. 

INPUT SPIKES MOST LIKELY DO NOT CONSTITUTE A LARGE FRACTION 
OF THE OBSERVED ACTIVITY 

The fraction of "input spikes" (drive) among all the spikes of the 
model is negligible at criticality (a = 1). This fraction, given a 



constant spike rate r, increases with tuning toward sub-criticality 
(a < 1), until all spikes are input spikes in the Poisson state 
(a = 0), and none arises from synaptic transmission. For exam- 
ple, in the driven, slightly sub-critical model (a = 0.99), only one 
in ~3600 spikes was an input spike. To illustrate this number, 
imagine a neuron that spikes with a rate of 1 Hz. This neuron 
fires spontaneously (i.e., an "input spike") only once an hour. 
This example is simplistic, because it assumes that the input is 
homogeneous, however, it illustrates well that the fraction of 
input spikes (from the external world, other brain structures, or 
of stochastic origin; see above) in the driven, sub-critical regime 
that reproduced the in vivo findings is extremely small compared 
to the overall level of activity. 

CONCEPTUAL CONSIDERATIONS ON THE ANALYSIS OF NEURAL 
AVALANCHES AND THE CRITICAL STATE 

While we have so far discussed how in vivo spike avalanches sug- 
gest a driven sub-critical regime of operation for mammalian ner- 
vous systems, several neglected but important conceptual issues 
with the analysis of neural avalanches surfaced in this study. These 
are discussed in the following. 

THE TERM "AVALANCHE" REFERS TO DIFFERENT ENTITIES IN SOC 
MODELS AND IN THE ANALYSIS OF NEURAL DATA 

Although it is rarely fully acknowledged, the term "avalanche" 
refers to different entities for activity in SOC models and for 
neural activity. In SOC models, an avalanche is a cascade of 
events that originates from a single input event (Bak et al., 1987). 
Subsequent avalanches are always separated by pauses (STS). In 
contrast, for neural activity, avalanches are defined using tem- 
poral binning (Beggs and Plenz, 2003), because neural activity 
lacks clear pauses that could naturally serve to define the begin- 
ning and end of an avalanche. Such avalanches can be defined 
on any spike time series, irrespective of its origin. Consequently, 
"binning-dependent avalanches" do not correspond to classical 
SOC avalanches. Although these two types of "avalanches" are 
different entities, it is customary to use the same term when 
referring to any one of them. In the present study, we analyzed 
the "binning-dependent" avalanches in both cases, in the driven 
models and in the in vivo activity. This justifies a comparison 
between model and in vivo activity, and was also necessary as 
binning-dependent avalanches are the de-facto standard in the 
analysis of neural systems, although previous studies frequently 
alluded to SOC avalanches. 

AVALANCHE DEFINITIONS IN HIGHLY PARALLEL RECORDINGS 

For neural activity, avalanches are commonly defined using tem- 
poral binning, and this definition relies on pauses. We can expect 
that physiologically relevant pauses (i.e., pauses of a few ms) van- 
ish in spike recordings, when activity of a large number of neurons 
is recorded in parallel. For example, if each neuron spikes with 
1 Hz, sampling only 100 neurons in parallel would frequently pro- 
duce pauses that are several milliseconds long. However, when 
sampling thousands or even millions of such neurons, pauses 
would probably be absent. Without pauses, neither the classical 
nor the binning-dependent avalanche definition is applicable, and 
consequently, alternative approaches to assess criticality have to 
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be established. Currently, these approaches threshold the activ- 
ity and thereby introduce pauses (e.g., Spasojevic et al., 1996; 
Papanikolaou et al., 2011; Poil et al., 2012). As an alternative 
approach, we propose to apply systematic subsampling. Both 
approaches allow using the binning-dependent avalanche defini- 
tions again. 

CAN WE DETERMINE A SPECIFIC CRITICAL EXPONENT FOR NEURAL 
DATA? 

Avalanche size distributions of critical branching processes have 
an exponent of r = 1.5 (Harris, 1963). Since branching processes 
have some resemblance with propagation of neural activity, it was 
hypothesized that neural avalanches should also show r = 1.5. 
Indeed, r = 1.5 has been observed (Beggs and Plenz, 2003, 2004; 
Stewart and Plenz, 2008; Hahn et al., 2010; Priesemann et al., 
2013), but only for specific bin sizes. For example, Beggs and 
Plenz showed in their seminal work that r = 1.5 holds for one 
specific bin size (4 ms), but when changing the bin size from 1 
to 16 ms, the exponent decreased from 2 to 1.2 (Beggs and Plenz, 
2003). Similarly, for the LFP avalanches shown here, x = 1.5 was 
observed only for a bin size of ~80ms, and with varying the 
bin size from 2.5 ms to ~2.5 s, the exponent changed from 3 
to 1 (Figure 10A) (Priesemann et al., 2013). Changes in x were 
also observed in the driven, subsampled SOC model (Figure 2C). 
Thus, drive and subsampling may underlie the variation of r 
in experiments as well. However, irrespective of its origin, it is 
an open question how to reconcile the variation of r in neu- 
ral data with the fixed x in critical systems. One proposal is to 
use a specific bin size for neural data, namely one average inter- 
event-interval (<IEI>) (Beggs and Plenz, 2003). However, there 
is no theoretical underpinning yet why this bin size should be 
preferred over others, and even for using this bin size, r was 
found to be 1.8 in spike avalanches in anesthetized cats (Hahn 
et al., 2010), instead of 1.5. Thus, in neural data, there is not 
a unique x, and therefore there is no specific critical exponent 
for neural activity, which would allow to link neural activity to a 
universality class. 

Since neural avalanche distributions change with the bin size 
(Beggs and Plenz, 2003; Priesemann et al, 2009, 2013; Benayoun 
et al., 2010; Hahn et al., 2010), we side with Benayoun et al., who 
"do not read any significance into the particular slope observed. 
[...] In our view, any good model of neural avalanches must 
reproduce the variability in the observed slope of the power law 
with temporal bin width." (Benayoun et al., 2010) Though we 
here did not observe power laws for the in vivo f(s), our model 
could reproduce the in vivo spike f(s) and their change with tem- 
poral binning. It could also reproduce the bin-size dependent 
changes of novel and established measures of avalanche dynam- 
ics (/(s = 1, bs), <s>, a*, DFA exponent). To the best of our 
knowledge, this is the first model that matched not only the 
avalanche properties for a single bin size, but also their changes 
with changing bin size. 

SUBSAMPLING EFFECTS IN THE ASSESSMENT OF CRITICALITY 

Subsampling is unavoidable in spike avalanche recordings in vivo, 
and is helpful when comparing neural activity to model activity 
(Priesemann et al., 2009). However, subsampling was also shown 



to complicate criticality analysis because it can distort avalanche 
measures (Priesemann et al., 2009, 2013; Ribeiro et al., 2010). 
To overcome this problem, we here developed avalanche mea- 
sures that are not distorted by subsampling. One example is the 
bin size dependence of the frequency of avalanches of size one 
(f{s = 1, bs)). This measure robustly showed power-law scal- 
ing in the driven SOC states, and exponential scaling in the 
Poisson state, independent of subsampling strategies (Figure 6). 
Therefore, we propose to use/(s =1, bs) as a robust measure for 
criticality analysis. 

Subsampling effects can appear very strong if one uses a fixed 
bin size, e.g., 1ms as in Ribeiro et al. (2014). We used instead 
a normalized bin size, which accounts for the problem that the 
population rate R changes with the number of sampled neurons. 
Using a normalized bin size diminished subsampling effects, and 
also allowed for a comparison to the in vivo recordings. 

FINITE SIZE EFFECTS IN CRITICALITY ASSESSMENT 

The finite size of the critical models limited the correlation 
lengths in space and time and thereby caused the cutoff in f(s) 
(Figures 2A,G). In analogy, the finite size is expected to also have 
caused - in the driven critical models - the cutoff at large bin size 
in the novel measures (f(s=l, bs), <s>). Since finite size effects 
decrease with increasing system size, and since the in vivo spikes 
were recorded in a far larger system than our model spikes, finite 
size effects are unlikely to account for the deviations from power 
law scaling found for the in vivo activity. 

In critical models, the finite size can change the value of a, 
for which the model is critical. For example, Eurich et al. (2002) 
showed for their model that the critical a depended on the model 
size L as a alt = 1 — L~ 0 5 . Thus, their finite size models with a — > 
1 were super-critical and showed peaks in their f(s). This was not 
the case for our critical models. Our models, in contrast, appeared 
to be slightly sub-critical at a = 1. This is probably due to the 
open boundary conditions we used in contrast to Eurich et al. 
Hence, since the finite size made our models at most sub-critical 
but not super-critical, there is no concern that the observed match 
of model and in vivo results at values of a < 1 is due to finite size 
effects. 

DIFFERENT TYPES OF CRITICAL PHASE TRANSITIONS EXIST 

To better understand criticality and potential deviations from it, 
it is also important to define which type of criticality one refers 
to. Critical phase transitions can occur for example for the tran- 
sitions from order to chaos (Bertschinger and Natschlager, 2004; 
Haldeman and Beggs, 2005; Boedecker et al., 2012; Lizier, 2013), 
from non-oscillatory to oscillatory regimes (Linkenkaer-Hansen 
et al., 2001; Poil et al., 2012), from replay to non-replay of spatio- 
temporal patterns (Scarpetta and de Candia, 2013), and from a 
regime with finite to one with potentially infinite avalanche sizes 
(Bak et al., 1987; Drossel and Schwabl, 1992; Olami et al, 1992; 
Eurich et al, 2002; Beggs and Plenz, 2003; Haldeman and Beggs, 
2005; Levina et al, 2007a,b, 2009), as known from branching 
processes (Harris, 1963). One study has found that the transitions 
to chaos and to potentially infinite avalanches coincide in their 
model (Haldeman and Beggs, 2005), but it is unclear whether this 
finding generalizes to other systems. We here want to emphasize 
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that our model showed a transition to potentially infinitely large 
avalanches. 

CONSEQUENCES FOR INFORMATION PROCESSING AND STABILITY OF 
BRAIN DYNAMICS 

After having discussed evidence from in vivo spike avalanche dis- 
tributions for a driven, sub-critical mode of operation, and after 
having clarified conceptual issues, we now turn to the question 
of what consequences these findings may have on information 
processing and dynamic stability in the mammalian brain. 

SUB-CRITICALITY, SUPER-CRITICALITY, AND STABILITY 

Criticality is characterized by a power-law distribution of its 
avalanche sizes. This indicates that avalanches of any size can 
occur; even close to infinite-size avalanches may occur, provided 
that the system is large enough to sustain them. Infinite-size 
avalanches do occur in the super-critical regime, and have been 
linked to epileptic seizures (Hsu et al., 2008; Meisel et al., 2012). 
Such infinite avalanches produce runaway activity, and could 
thereby impair normal brain activity. Therefore, it is unlikely that 
it would be good for a normally functioning brain to be super- 
critical. Sub -criticality, in contrast, never shows infinitely large 
avalanches, and thus offers a safer regime for brain operation. 
Thus, a slightly sub-critical regime allows the brain to avoid run- 
away activity, while still allowing moderate activity propagation, 
and maintaining most of the possible computational advantages 
that come with criticality (Haldeman and Beggs, 2005; Kinouchi 
and Copelli, 2006; Beggs, 2008; Shew et al, 2009; Shew and Plenz, 
2013). 

DRIVE AND INFORMATION PROCESSING 

There may be good reason why neural activity in vivo does 
not show a STS for its avalanches: When eliminating the STS, 
avalanches run in parallel, meet, and intermingle. Thereby, the 
rate of computations may be increased compared to the SOC 
state. In addition, the presence of multiple, potentially interact- 
ing avalanches, may enable collision-based computation, which is 
one fundamental way of information modification (Lizier, 2013). 
Thus, a driven state may increase the rate and capacity of neural 
information processing in vivo. 

CONCLUSIONS 

Our analysis of in vivo data indicated that the mammalian brain 
is not SOC because in vivo spiking activity differed fundamentally 
from activity expected for SOC. Instead, the mammalian brain 
apparently self-organizes to a slightly sub-critical regime with- 
out an STS. Mechanistically, such a driven, sub-critical regime 
shows a melange of avalanches, while SOC systems, in contrast, 
are characterized by temporally separated avalanches. Operating 
in a slightly sub-critical regime may prevent the brain from tip- 
ping over to super-criticality, which has been linked to epilepsy. 
Regarding computational capabilities, which have been reported 
to be optimal for SOC, a slightly sub-critical regime only deviates 
little from SOC and therefore its computational capabilities may 
still be close to optimal, while the non-zero drive in general may 
allow for a higher rate of information processing. Taken together, 
a driven, slightly sub-critical regime may strike a balance between 



optimal information processing and the need to avoid runaway 
activity. 

METHODS 

SELF-ORGANIZED CRITICAL MODEL 

The SOC neural network model we used here is the Bak-Tang- 
Wiesenfeld model (Bak et al, 1987), and modified versions of 
it. Translated to a neuroscience context, the model consisted of 
2500 non-leaky integrate and fire neurons. A neuron i spiked if its 
membrane voltage V,(t) reached a threshold 0: 

IfV,(t)>0, Vi(f+l) = Vj(t)-4 

© was set to © = 0 for convenience. Note that the choice of 
0 does not change the activity of the model at all. The model 
neurons were arranged on a 2D lattice, and each neuron was con- 
nected locally to its four next neighbors, i.e., the coupling strength 
otij = a for all four next neighbors of neuron i, and a;j = 0 else. 

Vi(t + 1) = Vi(t) + J2 a 'J- S( - f ~ + ff W 
i 

The time t was updated in ms (i.e., 1 ms effective synaptic delay). 
Tj denoted the spike times of neuron j, and H(t) was a func- 
tion which set a neuron above threshold with a certain Poisson 
rate h. h represented the "drive" in the context of SOC. Note 
that the neurons at the edges and corners of the grid had only 

3 and 2 neighbors, respectively. This model is equivalent to the 
well-known Bak- Tang- Wiesenfeld model (Bak et al., 1987) if 
h — > 0 and a = 1. In contrast, for a = 0, the model represented 
independent Poisson units which spiked with rate r = h. 

Subsampling (Priesemann et al, 2009) was applied to the 
model by sampling the activity of 100 randomly selected neu- 
rons only, and neglecting the activity of all other neurons. To 
simulate specific subsampling effects, the sampled neurons were 
not chosen randomly, but arranged in specific configurations (see 
Figure 6, right part). Here the sampled neurons were arranged to 
have very small or very large distances. For the small distances, 

4 x 4 or 8 x 8 neurons from a compact, central subset were sam- 
pled (Figure 6, red and pink), and for the large distances, 4x4 
or 8 x 8 neurons with distance 5 grid units between them were 
sampled (Figure 6, turquoise and beige). 

STOCHASTIC BRANCHING MODEL 

In addition to the SOC model, we also simulated a classical 
stochastic branching model. In this model, a branching process 
(Harris, 1963; Haldeman and Beggs, 2005) was mapped on a grid 
of neurons. An active neuron activated each of its k postsynap- 
tic neurons with probability p = a ■ l/k. As in the SOC model, 
this model was critical for a = 1 in the infinite size limit, and 
subcritical (supercritical) for a < 1 (a > 1). In contrast to the 
SOC model, here the postsynaptic neurons were assigned ran- 
domly at each step. The other parameters were analogous to the 
SOC model: The model had 2500 neurons with k = 4 connections 
each, and a and h were balances such that neurons spiked with 
r = 5 Hz (except if h — >■ 0). The open boundary conditions were 
implemented by defining pdiss = 0.001 as the probability that a 
neuron projected "outside of the grid," i.e., the probability that an 



Frontiers in Systems Neuroscience 



www.frontiersin.org 



June 2014 | Volume 8 | Article 108 | 13 



Priesemann et al. 



Spike avalanches in vivo 



activation of a postsynaptic neuron was not effective. Note that 
Pdiss > 0 makes the model slightly subcritical. Subsampling was 
implemented in the same manner as in the SOC model. Note 
however that spatial distances have no meaning in this model 
because of its random topology. Results for this model were qual- 
itatively similar to those of the SOC model. Therefore, we usually 
reported the results of the SOC model only. 

EXPERIMENTS 

We evaluated spikes from recordings in three different species, 
namely in rats, cats and monkeys. The rat experimental proto- 
cols were approved by the Institutional Animal Care and Use 
Committee of Rutgers University (Mizuseki et al., 2009). The 
cat experiments were performed in accordance with guidelines 
established by the Canadian Council for Animal Care (Blanche, 
2009). The monkey experiments were performed according to the 
German Law for the Protection of Experimental Animals, and 
were all approved by the Regierungsprasidium Darmstadt. The 
procedures also conformed to the regulations issued by the NIH 
and the Society for Neuroscience. 

The spike recordings from the rats and the cats came from 
the NSF-founded CRCNS data sharing website (Blanche, 2009; 
Mizuseki et al, 2009). In brief, in rats the spikes were recorded in 
CA1 of the right dorsal hippocampus during an open field task. 
We used the first data set of each animal (ec013.527, ec014.277, 
ec015.041, ec016.397), and from rat "ec014" we also used a sec- 
ond data set (ec014.333). The five datasets provided sorted spikes, 
i.e., {37, 77, 32, 58, 58} single units and {4, 8, 8, 8, 8} multi units, 
respectively. However, since the identity of a unit does not matter 
for the definition of neural avalanches (see below), the single- and 
multi-unit activity was combined to one set of spike times. More 
details on the experimental procedure and the datasets proper can 
be found on Mizuseki et al. (2009). 

For the spikes from the cat, neural data were recorded by Tim 
Blanche in the laboratory of Nicholas Swindale, University of 
British Columbia, and downloaded from the NSF-funded CRCNS 
Data Sharing website (Blanche, 2009). We used the data set pvc3, 
i.e., recordings in area 18 which contain 50 sorted single units 
(Blanche and Swindale, 2006). We used that part of the experi- 
ment in which no stimuli were presented, i.e., the spikes reflected 
spontaneous activity in the visual cortex of the anesthetized cat. 
Details on the experimental procedures and the data proper can 
be found in Blanche and Swindale (2006); Blanche (2009). 

In the monkey experiments, spikes were recorded simultane- 
ously from up to 16 single-ended micro -electrodes (0 = 80 fim) 
or tetrodes (0 = 96 /xm) in lateral prefrontal cortex of three 
trained macaque monkeys (Ml: 6kg $; M2: 12kg d; M3: 8kg 
9). The electrodes had impedances between 0.2 and 1.2 Mf2 at 
1 kHz, and were arranged in a square grid with inter electrode dis- 
tances of either 0.5 or 1.0 mm. The monkeys performed a visual 
short term memory task with on average 80% correct behav- 
ioral responses which required them to memorize a sample object 
and to compare a test stimulus presented after a delay of 3 s to 
memory content. The monkeys indicated via differential button 
press whether test and sample stimuli matched or not. Each trial 
consisted of a 1 s long baseline, 500-900 ms sample stimulus pre- 
sentation, a delay of 3 s and a response interval lasting throughout 



a 2 s test stimulus presentation. More details of the experimental 
procedure can be found in Pipa et al. (2009). In total, we ana- 
lyzed spike data from 1 1 experimental sessions comprising almost 
12.000 trials. In Ml and M2 we recorded four sessions each, and 
in M3 we recorded 3 sessions. 6 out of 1 1 sessions were recorded 
with tetrodes (2/4, 4/4, and 0/3 from Ml, M2, and M3, respec- 
tively). Spike sorting on the tetrode data was performed using 
a Bayesian optimal template matching approach as described in 
Franke (2011) (see Franke et al, 2010 for an earlier version) 
using the "Spyke Viewer" software (Propper and Obermayer, 
2013). On the single electrode data, spikes were sorted with a 
multi-dimensional PCA method (Smart Spike Sorter by Nan-Hui 
Chen). 

MEASURES 

Avalanches in SOC systems are cascades of spikes triggered by a 
single external spike (Bak et al., 1987). An avalanche can span 
the entire system, but can also affect just a few sites before it dies 
out. By definition, in SOC models subsequent avalanches are sep- 
arated by pauses that are much longer than the avalanches proper 
(STS) (Bak et al., 1987; Pruessner, 2012). This means that a new 
avalanche is only triggered after the previous one has long died 
out. In SOC systems, several avalanche characteristics, such as the 
distribution of sizes and durations, follow scaling laws, known 
from the framework of "renormalization theory" (Stanley, 1971, 
1999; Sethna et al., 2001; Dhar, 2006). In the following, we define 
the avalanche measures and describe the expected scaling laws for 
the SOC model and the critical stochastic branching model. 

The avalanche size s is the total number of spikes in an 
avalanche. The avalanche size distribution f(s) is its frequency 
of avalanche sizes, and p(s) refers to the respective probability 
distributions. f(s) follows a power law in SOC systems: 

f(s) ~ s- r 

r is the critical exponent and depends on the SOC model. For the 
SOC model we use here (a = 1 and h 0), x 1 (Bak et al., 
1987; Priesemann et al., 2009), and for the critical branching 
model r = 1.5 (Harris, 1963; Haldeman and Beggs, 2005). 

The definition of avalanche sizes in the driven models (h > 0) 
and in vivo relied on temporal binning (Beggs and Plenz, 2003), 
since these systems lacked STS. When applying temporal bins to a 
spike train, the avalanche size was defined as the total number of 
events in subsequent, non-empty time bins (Figure 1). Stating it 
differently, an avalanche is by definition the activity in a sequence 
of full bins, and is preceded and followed by an empty bin. With 
this definition, f(s) changed with the bin size (Figure 1). 

As stated above, f(s) changed with the bin size. To quan- 
tify the bin-size dependent changes of f(s), we used the mean 
avalanche size (<s>), and the measure /(s =1, bs), i.e., the bin 
size dependence of the frequency of avalanches of size s = 1 . 

A common measure to characterize neural avalanches is the 
branching parameter. In a branching process, the branching 
parameter a defines whether activity expands (a > 1) or dies out 
(er < 1) (Harris, 1963). Between these two regimes, at a = 1, the 
branching process is critical (Harris, 1963). In analogy, the a* 
was estimated from spike trains using temporal binning as follows 
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(Beggs and Plenz, 2003; Priesemann et al, 2009): a* is the num- 
ber of events in time bin f, divided by the number of events in 
time bin f,_i. The average over all a* (for which the number 
of events in f, _ j is not zero) is defined as the estimated branch- 
ing parameter a* (Figure 1) (Beggs and Plenz, 2003; Priesemann 
et al., 2009). Note that a* depends on the bin size, and may fail to 
provide the intended results (see Results and Discussion). 

Detrended fluctuations analysis (DFA) (Peng et al., 1994, 
1995; Kantelhardt et al., 2002) quantifies long-range correla- 
tions in a time-series, which also dominate SOC systems. We 
applied DFA to the time course of the summed population activ- 
ity. The summed population activity is the total number of 
spikes across all neurons at each sampling step. For the DFA, we 
used analysis window widths from 2 4 to 2 11 ms. Smaller win- 
dow widths could not be used because of the limited sampling 
resolution, and for windows larger than 2 s the power law scal- 
ing broke down, and this impeded the estimation of the DFA 
exponent /j. 

It sometimes is helpful to measure the bin size not in abso- 
lute time (e.g., milliseconds), but in "average inter event intervals" 
(<IEI>). The <IEI> is the inverse of the population rate _R, i.e., 
the rate of all units together, independent of their origin. In con- 
trast to the population rate _R, the rate of a single unit is denoted 
with r. 

LFP RECORDINGS IN HUMANS 

We evaluated LFP which were recorded with intracranial depth 
recordings in humans. We used the very same data and analysis 
methods as in Priesemann et al. (2013), and we used the results 
from all vigilance states combined, because we already showed 
that the differences with vigilance states were small (Priesemann 
et al., 2013). We analyzed data from five subjects [3 females (aged 
21, 23, and 27), two males; (aged 25 and 48)] with refractory 
partial epilepsy undergoing pre-surgical evaluation. The sub- 
jects were hospitalized between February 2005 and March 2007 
in the epilepsy unit at the Pitie-Salpetriere hospital in Paris. 
All patients gave their informed consent and procedures were 
approved by the local ethical committee (CCP). Each patient was 
continuously recorded during several days (duration range: 9-20 
days; mean duration: 16 days) with intracranial and scalp elec- 
trodes (Nicolet acquisition system, CA, US). Depth electrodes 
were composed of 4-10 cylindrical contacts (2.3-mm long, 1- 
mm in diameter, 10-mm apart center-to-center), mounted on 
a 1 mm wide flexible plastic probe. Pre and post implantation 
MRI scans were evaluated to anatomically locate each contact 
along the electrode trajectory. The placement of electrodes within 
each patient was determined solely by clinical criteria. Signals 
were digitized at 400 Hz. The five subjects were implanted with 
(44, 48, 50, 50, and 63) intracranial LFP recording sites. In total 
seven recording sites were excluded from the analysis due to 
artifacts and thus we used (44, 48, 45, 50, and 61) recordings 
sites for data evaluation. All LFP were low-pass filtered at 40 Hz 
(4th order butterworth, MATLAB) to reduce the impact of line 
noise. 

To analyze the neuronal avalanches for these LFP data in 
the same manner as the spike data, we extracted binary events 
from the LFP. These binary events represent phases of enhanced 



synaptic activity. To extract these events, we calculated the area 
under the positive deflection lobes between two zero crossings of 
the LFP (Figure 2 in Priesemann et al., 2013). As LFP- voltages 
reflect current flows via Ohm's law, this time integral, or area 
under the voltage curve, is proportional to the total amount 
of displaced charges and hence describes the departure from 
equilibrium (charge neutrality) quantitatively — in contrast to 
simple voltage peaks. To obtain binary events from the LFP, we 
applied a threshold to the area values under the LFP deflec- 
tion lobe. The threshold was selected such that each recording 
site in each interval of constant vigilance state had the same 
event rate r = 1 /4 Hz. In contrast to our first paper with these 
data (Priesemann et al., 2013), we here used only one value 
for r, and combined the results for all vigilance states from 
wakefulness to deep sleep, since neither r nor the different vig- 
ilance states affected the results qualitatively (Priesemann et al., 
2013). 

For the avalanche analysis in the humans, we used a bin size 
either in units of average inter event intervals ( <IEI> ) or in ms. 
The <IEI> is a function of the event rate r and the number 
of electrode contacts N, <IEI> = II {r ■ N) = 1/R. Since r was 
fixed and N did not vary much across patients, the following 
approximation holds: 1 <IEI> ~ 80 ms. 
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